Cladophora Microbiota, 2014, 16S Amplicon Analysis

#source('http://bioconductor.org/biocLite.R')
#biocLite('phyloseq')
library(phyloseq)
library(ggplot2)
library(plyr)
library(dplyr)
library(Rmisc)
library(DESeq2)
library(doParallel)
library(vegan)
library(grid)
library(gridExtra)
library(reshape2)

# Load biom file. 
biom <- import_biom("OTU_table.biom", "~/Dropbox/clado-manuscript/Nephele/PipelineResults_NMEPINZ20QK1/nephele_outputs/tree.tre", parseFunction=parse_taxonomy_greengenes)
biom

# Load and merge sample metadata with read data (biom). 
sam.data <- read.csv(file="sample.data.csv", row.names=1, header=TRUE)
sam.data$Date <- as.factor(sam.data$Date)
sam.data$DateSite <- paste(sam.data$Date, sam.data$Site)
sample_data(biom) <- sam.data
sample_data(biom)

## Sample Data:        [52 samples by 6 sample variables]:
##        TreatmentGroup  Site Date                       Description
## C178N1          Early North  178 Sample of day 178 at site North 1
## C178P1          Early Point  178 Sample of day 178 at site Point 1
## C185P2          Early Point  185 Sample of day 185 at site Point 2
## C206N2           Late North  206 Sample of day 206 at site North 2
## C206P1           Late Point  206 Sample of day 206 at site Point 1
## C206P2           Late Point  206 Sample of day 206 at site Point 2
## C214P1           Late Point  214 Sample of day 214 at site Point 1
## C214P2           Late Point  214 Sample of day 214 at site Point 2
## C214P3           Late Point  214 Sample of day 214 at site Point 3
## C214S1           Late South  214 Sample of day 214 at site South 1
## C214S2           Late South  214 Sample of day 214 at site South 2
## C214S3           Late South  214 Sample of day 214 at site South 3
## C185P1          Early Point  185 Sample of day 185 at site Point 1
## C185P3          Early Point  185 Sample of day 185 at site Point 3
## C199P3           Late Point  199 Sample of day 199 at site Point 3
## C199S2           Late South  199 Sample of day 199 at site South 2
## C199S3           Late South  199 Sample of day 199 at site South 3
## C206N1           Late North  206 Sample of day 206 at site North 1
## C178P2          Early Point  178 Sample of day 178 at site Point 2
## C199N3           Late North  199 Sample of day 199 at site North 3
## C206S1           Late South  206 Sample of day 206 at site South 1
## C214N3           Late North  214 Sample of day 214 at site North 3
## C199N2           Late North  199 Sample of day 199 at site North 2
## C206N3           Late North  206 Sample of day 206 at site North 3
## C206S2           Late South  206 Sample of day 206 at site South 2
## C199N1           Late North  199 Sample of day 199 at site North 1
## C199P1           Late Point  199 Sample of day 199 at site Point 1
## C199S1           Late South  199 Sample of day 199 at site South 1
## C214N1           Late North  214 Sample of day 214 at site North 1
## C172P1          Early Point  172 Sample of day 172 at site Point 1
## C199P2           Late Point  199 Sample of day 199 at site Point 2
## C172N3          Early North  172 Sample of day 172 at site North 3
## C172S3          Early South  172 Sample of day 172 at site South 3
## C178S2          Early South  178 Sample of day 178 at site South 2
## C178P3          Early Point  178 Sample of day 178 at site Point 3
## C178S3          Early South  178 Sample of day 178 at site South 3
## C172N1          Early North  172 Sample of day 172 at site North 1
## C172S1          Early South  172 Sample of day 172 at site South 1
## C178N3          Early North  178 Sample of day 178 at site North 3
## C185N2          Early North  185 Sample of day 185 at site North 2
## C185N3          Early North  185 Sample of day 185 at site North 3
## C185S3          Early South  185 Sample of day 185 at site South 3
## C214N2           Late North  214 Sample of day 214 at site North 2
## C172P2          Early Point  172 Sample of day 172 at site Point 2
## C185S2          Early South  185 Sample of day 185 at site South 2
## C172P3          Early Point  172 Sample of day 172 at site Point 3
## C185N1          Early North  185 Sample of day 185 at site North 1
## C172N2          Early North  172 Sample of day 172 at site North 2
## C178S1          Early South  178 Sample of day 178 at site South 1
## C185S1          Early South  185 Sample of day 185 at site South 1
## C172S2          Early South  172 Sample of day 172 at site South 2
## C178N2          Early North  178 Sample of day 178 at site North 2
##        SampleID.1  DateSite
## C178N1     C178N1 178 North
## C178P1     C178P1 178 Point
## C185P2     C185P2 185 Point
## C206N2     C206N2 206 North
## C206P1     C206P1 206 Point
## C206P2     C206P2 206 Point
## C214P1     C214P1 214 Point
## C214P2     C214P2 214 Point
## C214P3     C214P3 214 Point
## C214S1     C214S1 214 South
## C214S2     C214S2 214 South
## C214S3     C214S3 214 South
## C185P1     C185P1 185 Point
## C185P3     C185P3 185 Point
## C199P3     C199P3 199 Point
## C199S2     C199S2 199 South
## C199S3     C199S3 199 South
## C206N1     C206N1 206 North
## C178P2     C178P2 178 Point
## C199N3     C199N3 199 North
## C206S1     C206S1 206 South
## C214N3     C214N3 214 North
## C199N2     C199N2 199 North
## C206N3     C206N3 206 North
## C206S2     C206S2 206 South
## C199N1     C199N1 199 North
## C199P1     C199P1 199 Point
## C199S1     C199S1 199 South
## C214N1     C214N1 214 North
## C172P1     C172P1 172 Point
## C199P2     C199P2 199 Point
## C172N3     C172N3 172 North
## C172S3     C172S3 172 South
## C178S2     C178S2 178 South
## C178P3     C178P3 178 Point
## C178S3     C178S3 178 South
## C172N1     C172N1 172 North
## C172S1     C172S1 172 South
## C178N3     C178N3 178 North
## C185N2     C185N2 185 North
## C185N3     C185N3 185 North
## C185S3     C185S3 185 South
## C214N2     C214N2 214 North
## C172P2     C172P2 172 Point
## C185S2     C185S2 185 South
## C172P3     C172P3 172 Point
## C185N1     C185N1 185 North
## C172N2     C172N2 172 North
## C178S1     C178S1 178 South
## C185S1     C185S1 185 South
## C172S2     C172S2 172 South
## C178N2     C178N2 178 North

# Normalize by relative abundance. 
biom.relabund <- transform_sample_counts(biom, function(x) x / sum(x))

# Ordination plot, k = 3. 
ordNMDS.k3 <- ordinate(biom.relabund, method="NMDS", distance="bray", k=3)

## Run 0 stress 0.07869986 
## Run 1 stress 0.07869949 
## ... New best solution
## ... Procrustes: rmse 0.0002305481  max resid 0.001105219 
## ... Similar to previous best
## Run 2 stress 0.07869962 
## ... Procrustes: rmse 0.0002017883  max resid 0.001019135 
## ... Similar to previous best
## Run 3 stress 0.07869844 
## ... New best solution
## ... Procrustes: rmse 0.0004553477  max resid 0.001521332 
## ... Similar to previous best
## Run 4 stress 0.08334584 
## Run 5 stress 0.07874264 
## ... Procrustes: rmse 0.002511844  max resid 0.01380174 
## Run 6 stress 0.07869964 
## ... Procrustes: rmse 0.0005900723  max resid 0.002571066 
## ... Similar to previous best
## Run 7 stress 0.08334478 
## Run 8 stress 0.07872072 
## ... Procrustes: rmse 0.001339293  max resid 0.005271881 
## ... Similar to previous best
## Run 9 stress 0.07870018 
## ... Procrustes: rmse 0.0006557135  max resid 0.002997354 
## ... Similar to previous best
## Run 10 stress 0.07869808 
## ... New best solution
## ... Procrustes: rmse 0.0001766944  max resid 0.0005861921 
## ... Similar to previous best
## Run 11 stress 0.07869856 
## ... Procrustes: rmse 0.0002775179  max resid 0.001302053 
## ... Similar to previous best
## Run 12 stress 0.08337886 
## Run 13 stress 0.07869809 
## ... Procrustes: rmse 0.0001472413  max resid 0.0003891141 
## ... Similar to previous best
## Run 14 stress 0.07869902 
## ... Procrustes: rmse 0.0003348988  max resid 0.001535 
## ... Similar to previous best
## Run 15 stress 0.07869916 
## ... Procrustes: rmse 0.0003981774  max resid 0.001889146 
## ... Similar to previous best
## Run 16 stress 0.08334287 
## Run 17 stress 0.08804081 
## Run 18 stress 0.07869856 
## ... Procrustes: rmse 0.0002547032  max resid 0.001142475 
## ... Similar to previous best
## Run 19 stress 0.08695913 
## Run 20 stress 0.07906707 
## ... Procrustes: rmse 0.00900781  max resid 0.03814892 
## *** Solution reached

ord.k3 <- plot_ordination(biom.relabund, ordNMDS.k3, shape="Site", color = "Date") + geom_point(size=2)
ord.k3 + labs(title = "Cladophora, 2014") + theme_bw() + scale_colour_hue(h=c(300, 500))

# PERMANOVA. 
df = as(sample_data(biom), "data.frame")
d = phyloseq::distance(biom, "bray")
clado.adonis = adonis(d ~ Date*Site, df)
clado.adonis

## 
## Call:
## adonis(formula = d ~ Date * Site, data = df) 
## 
## Permutation: free
## Number of permutations: 999
## 
## Terms added sequentially (first to last)
## 
##           Df SumsOfSqs MeanSqs F.Model      R2 Pr(>F)    
## Date       5    3.1945 0.63890 17.4143 0.41338  0.001 ***
## Site       2    1.5561 0.77804 21.2068 0.20136  0.001 ***
## Date:Site 10    1.7298 0.17298  4.7148 0.22384  0.001 ***
## Residuals 34    1.2474 0.03669         0.16142           
## Total     51    7.7278                 1.00000           
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

biom.rich.est <- estimate_richness(biom, measures = NULL)
biom.rich.est$SampleID.1 <- row.names(biom.rich.est)
biom.rich.est <- merge(biom.rich.est, sam.data, by = "SampleID.1")
biom.rich.est$Date <- as.character(biom.rich.est$Date)
biom.rich.est$Date <- as.numeric(biom.rich.est$Date)
head(biom.rich.est)

##   SampleID.1 Observed    Chao1  se.chao1      ACE   se.ACE  Shannon
## 1     C172N1     5139 6770.782  99.67941 7316.921 48.50633 5.444561
## 2     C172N2     5236 7017.650 107.98761 7420.840 47.98329 5.649933
## 3     C172N3     6890 8499.605  90.45120 9053.801 49.75824 5.980043
## 4     C172P1     2936 4855.835 138.81016 5427.422 47.24285 4.851552
## 5     C172P2     4833 7043.680 130.90171 7619.611 51.59008 5.306940
## 6     C172P3     3750 5099.247  96.43281 5368.421 41.19882 4.824680
##     Simpson InvSimpson    Fisher TreatmentGroup  Site Date
## 1 0.9874631   79.76426  896.0810          Early North  172
## 2 0.9892749   93.23962 1001.0057          Early North  172
## 3 0.9899773   99.77399 1361.3171          Early North  172
## 4 0.9767986   43.10086  520.1978          Early Point  172
## 5 0.9803746   50.95444  904.5831          Early Point  172
## 6 0.9400545   16.68183  655.5148          Early Point  172
##                         Description  DateSite
## 1 Sample of day 172 at site North 1 172 North
## 2 Sample of day 172 at site North 2 172 North
## 3 Sample of day 172 at site North 3 172 North
## 4 Sample of day 172 at site Point 1 172 Point
## 5 Sample of day 172 at site Point 2 172 Point
## 6 Sample of day 172 at site Point 3 172 Point

# Plot observed richness. 
biom.rich.est.obs <- summarySE(biom.rich.est, measurevar="Observed", groupvars=c("Date","Site"))
p.obs <- ggplot(biom.rich.est.obs, aes(x=Date, y=Observed, color = Site)) + geom_point() +  geom_errorbar(aes(ymin=Observed-se, ymax=Observed+se)) + geom_line() + scale_colour_hue(h=c(400, 120))
p.obs + theme_bw() + theme(axis.text.x = element_text(size = 10, angle = 45, hjust=1),axis.text.y = element_text(size = 10))

# Plot Shanon index richness. 
biom.rich.est.sha <- summarySE(biom.rich.est, measurevar="Shannon", groupvars=c("Date","Site"))
p.sha <- ggplot(biom.rich.est.sha, aes(x=Date, y=Shannon, color = Site)) + geom_point() +  geom_errorbar(aes(ymin=Shannon-se, ymax=Shannon+se)) + geom_line() + scale_colour_hue(h=c(400, 120))
p.sha + theme_bw() + theme(axis.text.x = element_text(size = 10, angle = 45, hjust=1),axis.text.y = element_text(size = 10))

# Find top 30 genera and subset biom.relabund. 
sort.genera <- sort(tapply(taxa_sums(biom.relabund), tax_table(biom.relabund)[, "Genus"], sum), TRUE)
top.genera <- sort.genera[1:30]
top.genera.list <- names(top.genera)
biom.relabund.subset = subset_taxa(biom.relabund, Genus %in% top.genera.list)
biom.relabund.subset.taxa <- subset_taxa(biom.relabund.subset, Genus %in% as.factor(top.genera.list))
biom.relabund.subset.taxa

## phyloseq-class experiment-level object
## otu_table()   OTU Table:         [ 9093 taxa and 52 samples ]
## sample_data() Sample Data:       [ 52 samples by 6 sample variables ]
## tax_table()   Taxonomy Table:    [ 9093 taxa by 8 taxonomic ranks ]
## phy_tree()    Phylogenetic Tree: [ 9093 tips and 9089 internal nodes ]

relabund.top.genera <- psmelt(biom.relabund.subset.taxa)
relabund.top.genera.genus <- relabund.top.genera%>%
  group_by(Sample, Genus)%>%
  mutate(GenusAbundance = sum(Abundance))%>%
  distinct(Sample, GenusAbundance, TreatmentGroup, Site, Date, Family, Genus)
head(relabund.top.genera.genus)

## Source: local data frame [6 x 9]
## Groups: Sample, Genus [6]
## 
##   Sample GenusAbundance TreatmentGroup   Site   Date          Family
##    <chr>          <dbl>         <fctr> <fctr> <fctr>          <fctr>
## 1 C178S2     0.14728018          Early  South    178 Crenotrichaceae
## 2 C206N2     0.07541493           Late  North    206      Thermaceae
## 3 C199S1     0.12408425           Late  South    199 Crenotrichaceae
## 4 C185S1     0.12728636          Early  South    185 Crenotrichaceae
## 5 C199S3     0.10486769           Late  South    199 Crenotrichaceae
## 6 C172S1     0.10894857          Early  South    172   Cytophagaceae
## # ... with 3 more variables: Genus <fctr>, Sample <chr>, Genus <fctr>

# Summary of genus abundance of top 30 genera. 
relabund.top.genera.genus.est <- summarySE(relabund.top.genera.genus, measurevar="GenusAbundance", groupvars=c("Site","Date", "Genus"))
head(relabund.top.genera.genus.est)

##    Site Date         Genus N GenusAbundance           sd           se
## 1 North  172   Armatimonas 3   0.0082179013 0.0037601213 0.0021709070
## 2 North  172  Bdellovibrio 3   0.0023540936 0.0001963925 0.0001133873
## 3 North  172    Cellvibrio 3   0.0058399987 0.0066066072 0.0038143264
## 4 North  172          CM44 3   0.0088329518 0.0005298073 0.0003058844
## 5 North  172    Crenothrix 3   0.0008092387 0.0004823472 0.0002784833
## 6 North  172 Dechloromonas 3   0.0019270556 0.0026318679 0.0015195096
##             ci
## 1 0.0093406591
## 2 0.0004878661
## 3 0.0164117220
## 4 0.0013161143
## 5 0.0011982168
## 6 0.0065379223

relabund.top.genera.genus.est$Date <- as.character(relabund.top.genera.genus.est$Date)
relabund.top.genera.genus.est$Date <- as.numeric(relabund.top.genera.genus.est$Date)

# Plot summary of genus abundance of top 30 genera. 
p <- ggplot(relabund.top.genera.genus.est, aes(x=Date, y=GenusAbundance, color = Site)) + geom_point() +  geom_errorbar(aes(ymin=GenusAbundance-se, ymax=GenusAbundance+se)) +  geom_line() + facet_wrap(~Genus, ncol = 3, scales="free_y") + scale_colour_hue(h=c(400, 120))
p + theme_bw() + theme(axis.text.x = element_text(size = 10, angle = 45, hjust=1),axis.text.y = element_text(size = 10))

# Find methanotrophic bacteria by genus. 
methanolist <- read.table(file = "taxa-of-interest/methanos.txt")
methanolist <- as.vector(methanolist$V1)
biom.relabund.methanos <- subset_taxa(biom.relabund, Genus %in% as.factor(methanolist))
biom.relabund.methanos

## phyloseq-class experiment-level object
## otu_table()   OTU Table:         [ 567 taxa and 52 samples ]
## sample_data() Sample Data:       [ 52 samples by 6 sample variables ]
## tax_table()   Taxonomy Table:    [ 567 taxa by 8 taxonomic ranks ]
## phy_tree()    Phylogenetic Tree: [ 567 tips and 566 internal nodes ]

relabund.methanos <- psmelt(biom.relabund.methanos)
relabund.methanos.genus <- relabund.methanos%>%
  group_by(Sample, Genus)%>%
  mutate(GenusAbundance = sum(Abundance))%>%
  distinct(Sample, GenusAbundance, TreatmentGroup, Site, Date, Family, Genus)
relabund.methanos.genus.est <- summarySE(relabund.methanos.genus, measurevar="GenusAbundance", groupvars=c("Site","Date", "Genus"))
head(relabund.methanos.genus.est)

##    Site Date            Genus N GenusAbundance           sd           se
## 1 North  172       Crenothrix 3   8.092387e-04 4.823472e-04 2.784833e-04
## 2 North  172      Methylibium 3   6.600147e-04 3.800437e-04 2.194183e-04
## 3 North  172    Methylocaldum 3   7.717447e-05 5.988775e-05 3.457621e-05
## 4 North  172 Methylomicrobium 3   0.000000e+00 0.000000e+00 0.000000e+00
## 5 North  172     Methylomonas 3   1.097031e-04 2.736128e-06 1.579704e-06
## 6 North  172     Methylosinus 3   4.677898e-05 2.432451e-05 1.404376e-05
##             ci
## 1 1.198217e-03
## 2 9.440809e-04
## 3 1.487694e-04
## 4 0.000000e+00
## 5 6.796919e-06
## 6 6.042543e-05

relabund.methanos.genus.est$Date <- as.character(relabund.methanos.genus.est$Date)
relabund.methanos.genus.est$Date <- as.numeric(relabund.methanos.genus.est$Date)

# Plot summary of genus abundance of methanotrophic genera. 
p <- ggplot(relabund.methanos.genus.est, aes(x=Date, y=GenusAbundance, color = Site)) + geom_point() +  geom_errorbar(aes(ymin=GenusAbundance-se, ymax=GenusAbundance+se)) +  geom_line() + facet_wrap(~Genus, ncol = 3, scales="free_y") + scale_colour_hue(h=c(400, 120))
p + theme_bw() + theme(axis.text.x = element_text(size = 10, angle = 45, hjust=1),axis.text.y = element_text(size = 10))

# Find all genera
all.genera <- sort(get_taxa_unique(biom.relabund, "Genus"), decreasing=FALSE)
biom.relabund.all.genera <- subset_taxa(biom.relabund, Genus %in% as.factor(all.genera))
biom.relabund.all.genera <- psmelt(biom.relabund.all.genera)
biom.relabund.all.genera.genus <- biom.relabund.all.genera%>%
  group_by(Sample, Genus)%>%
  mutate(GenusAbundance = sum(Abundance))%>%
  distinct(Sample, GenusAbundance, TreatmentGroup, Site, Date, Family, Genus)

biom.relabund.all.genera.genus.est <- summarySE(biom.relabund.all.genera.genus, measurevar="GenusAbundance", groupvars=c("Site","Date", "Genus"))
head(biom.relabund.all.genera.genus.est)

##    Site Date           Genus N GenusAbundance           sd           se
## 1 North  172             A17 3   1.044934e-05 7.240333e-06 4.180208e-06
## 2 North  172   Achromobacter 3   1.205785e-06 2.088482e-06 1.205785e-06
## 3 North  172 Acidaminobacter 3   2.030206e-05 3.516421e-05 2.030206e-05
## 4 North  172      Acidocella 3   0.000000e+00 0.000000e+00 0.000000e+00
## 5 North  172      Acidovorax 3   4.023211e-05 3.508560e-05 2.025668e-05
## 6 North  172   Acinetobacter 3   1.395043e-04 1.803168e-04 1.041059e-04
##             ci
## 1 1.798598e-05
## 2 5.188076e-06
## 3 8.735273e-05
## 4 0.000000e+00
## 5 8.715747e-05
## 6 4.479317e-04

biom.relabund.all.genera.genus.est$Date <- as.character(biom.relabund.all.genera.genus.est$Date)
biom.relabund.all.genera.genus.est$Date <- as.numeric(biom.relabund.all.genera.genus.est$Date)
p <- ggplot(biom.relabund.all.genera.genus.est, aes(x=Date, y=GenusAbundance, color = Site)) + geom_point() +  geom_errorbar(aes(ymin=GenusAbundance-se, ymax=GenusAbundance+se)) +  geom_line() + facet_wrap(~Genus, ncol = 3, scales="free_y") + scale_colour_hue(h=c(400, 120))
p + theme_bw() + theme(axis.text.x = element_text(size = 10, angle = 45, hjust=1),axis.text.y = element_text(size = 10))

# We want to calculate the total relative abundance of each phylum. 
# "Melt" the phyloseq data into a dataframe and then take the top X most abundant phyla. 
biom.melt <- psmelt(biom.relabund)
biom.melt.sorted <- biom.melt %>%
  group_by(Sample,Phylum) %>%
  summarize(PhyAbund = sum(Abundance))%>%
  group_by(Phylum)%>%
  summarize(MeanPhyAbund = mean(PhyAbund))%>%
  arrange(-MeanPhyAbund)

# List of nPhyla top most abundant phyla. 
nPhyla =16
PhylumList <- biom.melt.sorted[1:nPhyla,1]
PhylumList <- PhylumList[is.na(PhylumList)==FALSE,]
PhylumList <- levels(droplevels(as.factor(PhylumList$Phylum)))
PhylumList

##  [1] "[Thermi]"         "Acidobacteria"    "Actinobacteria"  
##  [4] "Armatimonadetes"  "Bacteroidetes"    "Chlorobi"        
##  [7] "Chloroflexi"      "Cyanobacteria"    "Gemmatimonadetes"
## [10] "GN02"             "Planctomycetes"   "Proteobacteria"  
## [13] "SR1"              "TM7"              "Verrucomicrobia"

# Subset biom.melt for phyla. 
biom.subset <- subset_taxa(biom.relabund, Phylum %in% PhylumList)
biom.subset.melt <- psmelt(biom.subset)
biom.subset.melt.sorted = biom.subset.melt %>%
  group_by(Sample,Site,Date,Phylum) %>%
  summarize(PhyAbund = sum(Abundance))

# Summarize phylum abundances. 
biom.subset.melt.sorted.est <- summarySE(biom.subset.melt.sorted, measurevar="PhyAbund", groupvars=c("Site","Date", "Phylum"))
head(biom.subset.melt.sorted.est)

##    Site Date          Phylum N    PhyAbund           sd           se
## 1 North  172        [Thermi] 3 0.035290992 0.0143074443 0.0082604068
## 2 North  172   Acidobacteria 3 0.053632305 0.0245078281 0.0141496011
## 3 North  172  Actinobacteria 3 0.042155944 0.0102741604 0.0059317892
## 4 North  172 Armatimonadetes 3 0.008399230 0.0037857001 0.0021856750
## 5 North  172   Bacteroidetes 3 0.396831763 0.0276493911 0.0159633834
## 6 North  172        Chlorobi 3 0.000709788 0.0005910336 0.0003412334
##            ci
## 1 0.035541662
## 2 0.060880820
## 3 0.025522429
## 4 0.009404200
## 5 0.068684895
## 6 0.001468209

biom.subset.melt.sorted.est$Date <- as.character(biom.subset.melt.sorted.est$Date)
biom.subset.melt.sorted.est$Date <- as.numeric(biom.subset.melt.sorted.est$Date)

# Plot top phyla abundances. 
p <- ggplot(biom.subset.melt.sorted.est, aes(x=Date, y=PhyAbund, color = Site)) + geom_point() +  geom_errorbar(aes(ymin=PhyAbund-se, ymax=PhyAbund+se)) +  geom_line() + facet_wrap(~Phylum, ncol = 3, scales="free_y") + scale_colour_hue(h=c(400,120))
p + theme_bw() + theme(axis.text.x = element_text(size = 10, angle = 45, hjust=1),axis.text.y = element_text(size = 10)) + labs(x="Date",y="Relative Abundance") 

# Find top 30 classes and subset biom.relabund. 
sort.classes <- sort(tapply(taxa_sums(biom.relabund), tax_table(biom.relabund)[, "Class"], sum), TRUE)
top.classes <- sort.classes[1:30]
top.classes.list <- names(top.classes)
biom.relabund.subset = subset_taxa(biom.relabund, Class %in% top.classes.list)
biom.relabund.subset.taxa <- subset_taxa(biom.relabund.subset, Class %in% as.factor(top.classes.list))
biom.relabund.subset.taxa

## phyloseq-class experiment-level object
## otu_table()   OTU Table:         [ 38156 taxa and 52 samples ]
## sample_data() Sample Data:       [ 52 samples by 6 sample variables ]
## tax_table()   Taxonomy Table:    [ 38156 taxa by 8 taxonomic ranks ]
## phy_tree()    Phylogenetic Tree: [ 38156 tips and 38142 internal nodes ]

relabund.top.classes <- psmelt(biom.relabund.subset.taxa)
relabund.top.classes.class <- relabund.top.classes%>%
  group_by(Sample, Class)%>%
  mutate(ClassAbundance = sum(Abundance))%>%
  distinct(Sample, ClassAbundance, TreatmentGroup, Site, Date, Family, Class)
relabund.top.classes.class.print <- summarySE(relabund.top.classes.class, measurevar="ClassAbundance", groupvars=c("Class"))
relabund.top.classes.class.print[,c("Class","N","ClassAbundance")]

##                    Class    N ClassAbundance
## 1   [Chloracidobacteria]  156   0.0242780076
## 2         [Pedosphaerae]  312   0.0008576936
## 3          [Saprospirae]  156   0.2600740284
## 4         Acidimicrobiia  572   0.0142314956
## 5        Acidobacteria-6  156   0.0087173757
## 6         Actinobacteria 1092   0.0008919558
## 7    Alphaproteobacteria 1196   0.1189849423
## 8          Armatimonadia   52   0.0043355902
## 9            Bacteroidia  468   0.0016047794
## 10                 BD1-5   52   0.0036947445
## 11    Betaproteobacteria  728   0.1394166181
## 12           Chloroplast  156   0.0387461867
## 13            Clostridia  624   0.0008200079
## 14            Cytophagia  260   0.0479689304
## 15            Deinococci  156   0.0590344529
## 16   Deltaproteobacteria 1352   0.0199104462
## 17        Flavobacteriia  208   0.0494758417
## 18   Gammaproteobacteria 1664   0.0668946728
## 19      Gemmatimonadetes  156   0.0019508196
## 20                 OM190   52   0.0040714493
## 21                 OPB56   52   0.0007609416
## 22         Phycisphaerae  104   0.0016755586
## 23        Planctomycetia  260   0.0069005428
## 24          Solibacteres  208   0.0118832190
## 25      Sphingobacteriia  156   0.0073552263
## 26 Synechococcophycideae  260   0.0354155284
## 27                  TA18   52   0.0015786666
## 28                 TM7-1   52   0.0018488588
## 29      Verrucomicrobiae   52   0.0026770285
## 30                   ZB2   52   0.0007206090

# Summary of class abundance of top 30 classes. 
relabund.top.classes.class.est <- summarySE(relabund.top.classes.class, measurevar="ClassAbundance", groupvars=c("Site","Date", "Class"))
head(relabund.top.classes.class.est)

##    Site Date                Class  N ClassAbundance           sd
## 1 North  172 [Chloracidobacteria]  9    0.039306658 0.0176302507
## 2 North  172       [Pedosphaerae] 18    0.001325860 0.0001686127
## 3 North  172        [Saprospirae]  9    0.191751831 0.0190289005
## 4 North  172       Acidimicrobiia 33    0.038384578 0.0085857752
## 5 North  172      Acidobacteria-6  9    0.003818652 0.0006679931
## 6 North  172       Actinobacteria 63    0.003232369 0.0001236175
##             se           ci
## 1 5.876750e-03 1.355181e-02
## 2 3.974240e-05 8.384913e-05
## 3 6.342967e-03 1.462691e-02
## 4 1.494592e-03 3.044383e-03
## 5 2.226644e-04 5.134650e-04
## 6 1.557434e-05 3.113266e-05

relabund.top.classes.class.est$Date <- as.character(relabund.top.classes.class.est$Date)
relabund.top.classes.class.est$Date <- as.numeric(relabund.top.classes.class.est$Date)

# Plot summary of class abundance of top 30 classes. 
p <- ggplot(relabund.top.classes.class.est, aes(x=Date, y=ClassAbundance, color = Site)) + geom_point() +  geom_errorbar(aes(ymin=ClassAbundance-se, ymax=ClassAbundance+se)) +  geom_line() + facet_wrap(~Class, ncol = 3, scales="free_y") + scale_colour_hue(h=c(400, 120))
p + theme_bw() + theme(axis.text.x = element_text(size = 10, angle = 45, hjust=1),axis.text.y = element_text(size = 10))

# Find all classes
all.classes <- sort(get_taxa_unique(biom.relabund, "Class"), decreasing=FALSE)
biom.relabund.all.classes <- subset_taxa(biom.relabund, Class %in% as.factor(all.classes))
biom.relabund.all.classes <- psmelt(biom.relabund.all.classes)
biom.relabund.all.classes.class <- biom.relabund.all.classes%>%
  group_by(Sample, Class)%>%
  mutate(ClassAbundance = sum(Abundance))%>%
  distinct(Sample, ClassAbundance, TreatmentGroup, Site, Date, Family, Class)

biom.relabund.all.classes.class.est <- summarySE(biom.relabund.all.classes.class, measurevar="ClassAbundance", groupvars=c("Site","Date", "Class"))
head(biom.relabund.all.classes.class.est)

##    Site Date                Class N ClassAbundance           sd
## 1 North  172       [Brachyspirae] 6   0.000000e+00 0.000000e+00
## 2 North  172       [Brevinematae] 3   0.000000e+00 0.000000e+00
## 3 North  172 [Chloracidobacteria] 9   3.930666e-02 1.763025e-02
## 4 North  172        [Cloacamonae] 6   0.000000e+00 0.000000e+00
## 5 North  172     [Fimbriimonadia] 3   2.686822e-05 1.191337e-05
## 6 North  172      [Lentisphaeria] 6   1.561697e-06 2.419371e-06
##             se           ci
## 1 0.000000e+00 0.000000e+00
## 2 0.000000e+00 0.000000e+00
## 3 5.876750e-03 1.355181e-02
## 4 0.000000e+00 0.000000e+00
## 5 6.878186e-06 2.959445e-05
## 6 9.877040e-07 2.538974e-06

biom.relabund.all.classes.class.est$Date <- as.character(biom.relabund.all.classes.class.est$Date)
biom.relabund.all.classes.class.est$Date <- as.numeric(biom.relabund.all.classes.class.est$Date)

p <- ggplot(biom.relabund.all.classes.class.est, aes(x=Date, y=ClassAbundance, color = Site)) + geom_point() +  geom_errorbar(aes(ymin=ClassAbundance-se, ymax=ClassAbundance+se)) +  geom_line() + facet_wrap(~Class, ncol = 3, scales="free_y") + scale_colour_hue(h=c(400, 120))
p + theme_bw() + theme(axis.text.x = element_text(size = 10, angle = 45, hjust=1),axis.text.y = element_text(size = 10))

# Data downloaded from 
data <- read.csv("north_temperate_lakes_lter__daily_water_temperature_-_lake_mendota.csv", header=T)
head(data)

##   sampledate year4 month daynum depth  wtemp flag_wtemp
## 1 2006-06-27  2006     6    178   0.0 22.205          N
## 2 2006-06-27  2006     6    178   0.5 22.202          N
## 3 2006-06-27  2006     6    178   1.0 22.244          N
## 4 2006-06-27  2006     6    178   1.5 22.219          N
## 5 2006-06-27  2006     6    178   2.0 22.259          N
## 6 2006-06-27  2006     6    178   2.5 22.214          N

nolegend <- theme(legend.position="none")
Jdays <- c(172,178,185,199,206,214)

# Subset all wtemp at or below 1.5 meters (sample depth up to shore). 
clado.depth <- subset(data, depth<=1.5, select = year4:wtemp)

# All years (2006-2016), full year. 
p <- qplot(clado.depth$daynum, y = clado.depth$wtemp)
p <- p + geom_point(aes(colour = clado.depth$year4), size=1, alpha = 0.8) + 
  scale_y_continuous(limits = c(4,30))+
  labs(title = "Water Temperature <1.5m Depth, Lake Mendota") + 
  xlab("Julian Day") + ylab("Temperature °C")+ 
  geom_vline(xintercept=Jdays, colour="darkgreen", linetype = "longdash") +
  guides(color=guide_legend(title="Year")) + 
  scale_colour_gradient(low = "darkred")
p+ theme_bw() + theme(axis.text.x = element_text(size = 10, angle = 45, hjust=1),axis.text.y = element_text(size = 10))

# All years (2006-2016), sample dates ± 25 days
clado.depth.zoom <- subset(clado.depth, daynum>=150 & daynum<=250, select = year4:wtemp)
p <- qplot(clado.depth.zoom$daynum, y = clado.depth.zoom$wtemp) + 
  scale_y_continuous(limits = c(10,30))+
  geom_point(aes(colour = clado.depth.zoom$year4), size=1, alpha = 0.8) + 
  labs(title = "Water Temperature <1.5m Depth, Lake Mendota", x="Julian Day", y="Temperature °C") + 
  geom_vline(xintercept=Jdays, colour="darkgreen", linetype = "longdash") + 
  guides(color=guide_legend(title="Year")) + scale_colour_gradient(low = "darkred")
p+ theme_bw() + theme(axis.text.x = element_text(size = 10, angle = 45, hjust=1),axis.text.y = element_text(size = 10))

# Smarter way to make grid of plots. 
p <- ggplot(clado.depth.zoom, aes(daynum, wtemp))
p <- p + geom_point(size = 0.5) + 
  ylim(15,30) + 
  facet_wrap(~year4, ncol = 3) + 
  geom_vline(xintercept=Jdays, colour="grey30", linetype = "longdash") +
  guides(color=guide_legend(title="Year")) +
  labs(title = "Water Temperature <1.5m Depth, Lake Mendota", x="Julian Day", y="Temperature °C")
p + theme_bw() + theme(axis.text.x = element_text(size = 10, angle = 45, hjust=1),axis.text.y = element_text(size = 10))

# Find genera of interest and subset biom.relabund. 
genofint <- read.table(file = "taxa-of-interest/genera-of-interest.txt")
genofint.list <- as.vector(genofint$V1)
biom.relabund.subset.genofint = subset_taxa(biom.relabund, Genus %in% genofint.list)
biom.relabund.subset.genofint.taxa <- subset_taxa(biom.relabund.subset.genofint, Genus %in% as.factor(genofint.list))
biom.relabund.subset.genofint.taxa

## phyloseq-class experiment-level object
## otu_table()   OTU Table:         [ 6620 taxa and 52 samples ]
## sample_data() Sample Data:       [ 52 samples by 6 sample variables ]
## tax_table()   Taxonomy Table:    [ 6620 taxa by 8 taxonomic ranks ]
## phy_tree()    Phylogenetic Tree: [ 6620 tips and 6618 internal nodes ]

relabund.genofint <- psmelt(biom.relabund.subset.genofint.taxa)
relabund.genofint.genus <- relabund.genofint%>%
  group_by(Sample, Genus)%>%
  mutate(GenusAbundance = sum(Abundance))%>%
  distinct(Sample, GenusAbundance, TreatmentGroup, Site, Date, Family, Genus)
head(relabund.genofint.genus)

## Source: local data frame [6 x 9]
## Groups: Sample, Genus [6]
## 
##   Sample GenusAbundance TreatmentGroup   Site   Date          Family
##    <chr>          <dbl>         <fctr> <fctr> <fctr>          <fctr>
## 1 C178S2     0.14728018          Early  South    178 Crenotrichaceae
## 2 C206N2     0.07541493           Late  North    206      Thermaceae
## 3 C199S1     0.12408425           Late  South    199 Crenotrichaceae
## 4 C185S1     0.12728636          Early  South    185 Crenotrichaceae
## 5 C199S3     0.10486769           Late  South    199 Crenotrichaceae
## 6 C172S1     0.10894857          Early  South    172   Cytophagaceae
## # ... with 3 more variables: Genus <fctr>, Sample <chr>, Genus <fctr>

# Summary of genus abundance of genera of interest. 
relabund.genofint.genus.est <- summarySE(relabund.genofint.genus, measurevar="GenusAbundance", groupvars=c("Site","Date", "Genus"))
head(relabund.genofint.genus.est)

##    Site Date          Genus N GenusAbundance           sd           se
## 1 North  172    Armatimonas 3   0.0082179013 0.0037601213 0.0021709070
## 2 North  172   Bdellovibrio 3   0.0023540936 0.0001963925 0.0001133873
## 3 North  172     Crenothrix 3   0.0008092387 0.0004823472 0.0002784833
## 4 North  172 Flavobacterium 3   0.0853337233 0.0369814367 0.0213512424
## 5 North  172 Flectobacillus 3   0.0109642505 0.0028549117 0.0016482840
## 6 North  172     Fluviicola 3   0.0032900333 0.0012430598 0.0007176809
##             ci
## 1 0.0093406591
## 2 0.0004878661
## 3 0.0011982168
## 4 0.0918669814
## 5 0.0070919938
## 6 0.0030879318

relabund.genofint.genus.est$Date <- as.character(relabund.genofint.genus.est$Date)
relabund.genofint.genus.est$Date <- as.numeric(relabund.genofint.genus.est$Date)

# Plot summary of genus abundance of genera of interest. 
p <- ggplot(relabund.genofint.genus.est, aes(x=Date, y=GenusAbundance, color = Site)) + geom_point() +  geom_errorbar(aes(ymin=GenusAbundance-se, ymax=GenusAbundance+se)) +  geom_line() + facet_wrap(~Genus, ncol = 2, scales="free_y") + scale_colour_hue(h=c(400, 120))
p + theme_bw() + theme(axis.text.x = element_text(size = 10, angle = 45, hjust=1),axis.text.y = element_text(size = 10))

# Find classes_zulk of interest and subset biom.relabund. 
classzulk <- read.csv(file = "taxa-of-interest/classes-zulkifly.txt", header = T)
classzulk.list <- as.vector(classzulk$Class)
biom.relabund.subset.classzulk = subset_taxa(biom.relabund, Class %in% classzulk.list)
biom.relabund.subset.classzulk.taxa <- subset_taxa(biom.relabund.subset.classzulk, Class %in% as.factor(classzulk.list))
biom.relabund.subset.classzulk.taxa

## phyloseq-class experiment-level object
## otu_table()   OTU Table:         [ 24232 taxa and 52 samples ]
## sample_data() Sample Data:       [ 52 samples by 6 sample variables ]
## tax_table()   Taxonomy Table:    [ 24232 taxa by 8 taxonomic ranks ]
## phy_tree()    Phylogenetic Tree: [ 24232 tips and 24220 internal nodes ]

relabund.classzulk <- psmelt(biom.relabund.subset.classzulk.taxa)
relabund.classzulk.class <- relabund.classzulk%>%
  group_by(Sample, Class)%>%
  mutate(ClassAbundance = sum(Abundance))%>%
  distinct(Sample, ClassAbundance, TreatmentGroup, Site, Date, Family, Class)
head(relabund.classzulk.class)

## Source: local data frame [6 x 9]
## Groups: Sample, Class [6]
## 
##   Sample ClassAbundance TreatmentGroup   Site   Date          Family
##    <chr>          <dbl>         <fctr> <fctr> <fctr>          <fctr>
## 1 C206N3      0.2638815           Late  North    206  Aeromonadaceae
## 2 C178S2      0.1751127          Early  South    178 Crenotrichaceae
## 3 C206N2      0.1210454           Late  North    206      Thermaceae
## 4 C199S1      0.1578279           Late  South    199 Crenotrichaceae
## 5 C185S1      0.1641425          Early  South    185 Crenotrichaceae
## 6 C185N2      0.1079530          Early  North    185    Trueperaceae
## # ... with 3 more variables: Class <fctr>, Sample <chr>, Class <fctr>

# Summary of class abundance of classes_zulk of interest. 
relabund.classzulk.class.est <- summarySE(relabund.classzulk.class, measurevar= "ClassAbundance", groupvars=c("Site","Date","Class"))
head(relabund.classzulk.class.est)

##    Site Date               Class  N ClassAbundance          sd
## 1 North  172 Alphaproteobacteria 69    0.164994059 0.024535262
## 2 North  172         Bacteroidia 27    0.003037096 0.003668845
## 3 North  172  Betaproteobacteria 42    0.134184504 0.012226685
## 4 North  172          Deinococci  9    0.035290992 0.012390610
## 5 North  172 Deltaproteobacteria 78    0.015869340 0.001210423
## 6 North  172      Flavobacteriia 12    0.100069095 0.032296939
##             se           ci
## 1 0.0029536985 0.0058940135
## 2 0.0007060695 0.0014513466
## 3 0.0018866184 0.0038101032
## 4 0.0041302034 0.0095242661
## 5 0.0001370534 0.0002729082
## 6 0.0093233233 0.0205204961

relabund.classzulk.class.est$Date <- as.character(relabund.classzulk.class.est$Date)
relabund.classzulk.class.est$Date <- as.numeric(relabund.classzulk.class.est$Date)

# Plot MERGED summary of class abundance of classes_zulk of interest and Zulkifly et al. class abundances. 
relabund.classzulk.class.est.merge <- merge(relabund.classzulk.class.est, classzulk, by = "Class")
relabund.classzulk.class.est.merge

##                   Class  Site Date  N ClassAbundance           sd
## 1   Alphaproteobacteria North  172 69   1.649941e-01 2.453526e-02
## 2   Alphaproteobacteria Point  206 46   9.164773e-02 6.169728e-03
## 3   Alphaproteobacteria Point  199 69   9.852536e-02 1.085808e-02
## 4   Alphaproteobacteria Point  185 69   9.225076e-02 8.204794e-03
## 5   Alphaproteobacteria Point  178 69   1.041032e-01 5.367296e-03
## 6   Alphaproteobacteria South  214 69   1.278667e-01 7.796891e-03
## 7   Alphaproteobacteria Point  172 69   9.560667e-02 7.851879e-03
## 8   Alphaproteobacteria North  214 69   1.116651e-01 7.125619e-03
## 9   Alphaproteobacteria North  206 69   1.108972e-01 2.441773e-02
## 10  Alphaproteobacteria South  185 69   9.285924e-02 4.715822e-03
## 11  Alphaproteobacteria North  199 69   1.254696e-01 1.205607e-02
## 12  Alphaproteobacteria North  185 69   1.314984e-01 2.635229e-03
## 13  Alphaproteobacteria North  178 69   1.404780e-01 1.300084e-02
## 14  Alphaproteobacteria Point  214 69   9.559706e-02 8.922345e-04
## 15  Alphaproteobacteria South  206 46   1.638028e-01 6.242964e-03
## 16  Alphaproteobacteria South  199 69   1.003508e-01 5.144320e-03
## 17  Alphaproteobacteria South  178 69   1.389276e-01 1.247185e-02
## 18  Alphaproteobacteria South  172 69   1.610157e-01 3.121611e-03
## 19          Bacteroidia North  172 27   3.037096e-03 3.668845e-03
## 20          Bacteroidia Point  206 18   1.090664e-03 2.244005e-04
## 21          Bacteroidia Point  178 27   1.170437e-03 1.944241e-04
## 22          Bacteroidia South  214 27   4.882563e-03 2.665430e-04
## 23          Bacteroidia North  206 27   1.440295e-03 2.851560e-04
## 24          Bacteroidia South  185 27   5.683288e-04 1.061741e-04
## 25          Bacteroidia Point  199 27   7.969615e-04 1.290053e-04
## 26          Bacteroidia Point  185 27   4.259956e-03 9.481523e-04
## 27          Bacteroidia North  178 27   2.833241e-04 1.567245e-04
## 28          Bacteroidia Point  214 27   1.112868e-03 1.704668e-04
## 29          Bacteroidia Point  172 27   3.064871e-04 8.381442e-05
## 30          Bacteroidia North  214 27   1.999671e-03 2.125027e-04
## 31          Bacteroidia South  199 27   2.991933e-03 1.633787e-04
## 32          Bacteroidia North  199 27   1.577465e-03 1.583485e-04
## 33          Bacteroidia North  185 27   2.546528e-04 1.595386e-04
## 34          Bacteroidia South  172 27   6.550725e-04 3.207826e-04
## 35          Bacteroidia South  206 18   1.105381e-03 1.471926e-04
## 36          Bacteroidia South  178 27   1.015037e-03 2.778247e-04
## 37   Betaproteobacteria North  172 42   1.341845e-01 1.222668e-02
## 38   Betaproteobacteria Point  178 42   1.814319e-01 1.532576e-02
## 39   Betaproteobacteria North  206 42   1.041486e-01 2.104803e-02
## 40   Betaproteobacteria Point  206 28   9.966041e-02 4.315876e-03
## 41   Betaproteobacteria Point  199 42   1.091130e-01 9.293773e-03
## 42   Betaproteobacteria North  178 42   1.637763e-01 3.320126e-02
## 43   Betaproteobacteria Point  214 42   8.081258e-02 3.605772e-03
## 44   Betaproteobacteria South  214 42   1.420981e-01 3.528373e-03
## 45   Betaproteobacteria North  214 42   1.109697e-01 4.704533e-03
## 46   Betaproteobacteria South  199 42   1.522503e-01 6.007583e-03
## 47   Betaproteobacteria South  185 42   1.690466e-01 4.865455e-03
## 48   Betaproteobacteria Point  185 42   1.299196e-01 1.267453e-02
## 49   Betaproteobacteria North  185 42   1.585639e-01 1.197191e-02
## 50   Betaproteobacteria South  172 42   1.628000e-01 2.230144e-02
## 51   Betaproteobacteria Point  172 42   1.299117e-01 1.713028e-02
## 52   Betaproteobacteria South  206 28   1.186570e-01 7.628531e-03
## 53   Betaproteobacteria North  199 42   1.557116e-01 8.083082e-03
## 54   Betaproteobacteria South  178 42   1.862712e-01 2.762646e-02
## 55           Deinococci North  172  9   3.529099e-02 1.239061e-02
## 56           Deinococci Point  206  6   3.159299e-02 8.097845e-04
## 57           Deinococci North  178  9   6.884618e-02 1.342334e-02
## 58           Deinococci Point  214  9   5.218391e-02 8.947982e-04
## 59           Deinococci North  214  9   8.296689e-02 1.548303e-03
## 60           Deinococci Point  178  9   7.742805e-02 2.268726e-03
## 61           Deinococci Point  199  9   7.389939e-02 4.483164e-03
## 62           Deinococci Point  185  9   6.721340e-02 3.932361e-03
## 63           Deinococci North  185  9   9.413589e-02 1.207945e-02
## 64           Deinococci North  206  9   9.772030e-02 1.895029e-02
## 65           Deinococci South  199  9   4.399603e-02 5.515838e-03
## 66           Deinococci South  214  9   4.089510e-02 1.914282e-03
## 67           Deinococci Point  172  9   4.024737e-02 1.769149e-02
## 68           Deinococci South  172  9   2.589008e-02 1.067506e-02
## 69           Deinococci South  185  9   4.119949e-02 8.925328e-04
## 70           Deinococci North  199  9   8.387299e-02 1.078240e-02
## 71           Deinococci South  206  6   5.233920e-02 1.207813e-04
## 72           Deinococci South  178  9   4.152299e-02 1.185636e-02
## 73  Deltaproteobacteria North  178 78   8.647701e-03 1.174233e-03
## 74  Deltaproteobacteria North  185 78   8.745398e-03 1.265861e-03
## 75  Deltaproteobacteria Point  214 78   1.603615e-02 9.944642e-04
## 76  Deltaproteobacteria North  172 78   1.586934e-02 1.210423e-03
## 77  Deltaproteobacteria Point  206 52   2.362662e-02 1.830317e-03
## 78  Deltaproteobacteria South  172 78   2.068596e-02 2.030225e-03
## 79  Deltaproteobacteria Point  199 78   2.076148e-02 1.425853e-03
## 80  Deltaproteobacteria Point  185 78   1.215409e-02 1.129379e-03
## 81  Deltaproteobacteria North  199 78   2.700271e-02 4.455407e-03
## 82  Deltaproteobacteria South  178 78   1.178370e-02 1.839846e-03
## 83  Deltaproteobacteria Point  172 78   1.310480e-02 2.840589e-03
## 84  Deltaproteobacteria North  214 78   2.415750e-02 1.255331e-03
## 85  Deltaproteobacteria Point  178 78   9.133387e-03 1.666179e-03
## 86  Deltaproteobacteria South  214 78   3.072906e-02 1.439476e-03
## 87  Deltaproteobacteria North  206 78   3.789912e-02 1.649665e-02
## 88  Deltaproteobacteria South  185 78   1.791409e-02 8.036882e-04
## 89  Deltaproteobacteria South  206 52   3.416406e-02 7.891517e-04
## 90  Deltaproteobacteria South  199 78   3.196279e-02 2.115607e-03
## 91       Flavobacteriia North  178 12   3.939026e-02 3.074583e-03
## 92       Flavobacteriia Point  214 12   4.980656e-02 2.572204e-03
## 93       Flavobacteriia North  185 12   3.812167e-02 7.076609e-03
## 94       Flavobacteriia South  172 12   7.149890e-02 1.402093e-02
## 95       Flavobacteriia North  199 12   6.037744e-02 4.392248e-03
## 96       Flavobacteriia South  178 12   2.352850e-02 3.886989e-03
## 97       Flavobacteriia Point  199 12   4.979223e-02 3.629757e-03
## 98       Flavobacteriia North  172 12   1.000691e-01 3.229694e-02
## 99       Flavobacteriia Point  178 12   6.954441e-02 2.698822e-03
## 100      Flavobacteriia Point  206  8   6.637744e-02 5.019688e-03
## 101      Flavobacteriia Point  185 12   3.386490e-02 4.224023e-03
## 102      Flavobacteriia North  206 12   2.977471e-02 4.835590e-03
## 103      Flavobacteriia South  214 12   4.044778e-02 2.668235e-03
## 104      Flavobacteriia Point  172 12   7.058525e-02 4.880392e-03
## 105      Flavobacteriia North  214 12   3.770521e-02 3.077634e-03
## 106      Flavobacteriia South  185 12   2.764441e-02 4.584770e-04
## 107      Flavobacteriia South  199 12   4.273863e-02 1.894567e-03
## 108      Flavobacteriia South  206  8   4.265953e-02 1.466844e-03
## 109 Gammaproteobacteria North  178 96   3.280080e-02 3.213318e-03
## 110 Gammaproteobacteria North  185 96   3.261424e-02 4.501205e-04
## 111 Gammaproteobacteria North  199 96   2.336821e-02 3.449423e-04
## 112 Gammaproteobacteria South  172 96   2.336364e-02 1.509224e-03
## 113 Gammaproteobacteria Point  214 96   3.444340e-02 2.229327e-03
## 114 Gammaproteobacteria North  172 96   4.623499e-02 2.368133e-02
## 115 Gammaproteobacteria South  178 96   1.488996e-01 2.136030e-02
## 116 Gammaproteobacteria Point  199 96   3.140871e-02 2.329793e-03
## 117 Gammaproteobacteria Point  178 96   6.566709e-02 6.560287e-03
## 118 Gammaproteobacteria Point  206 64   4.262246e-02 1.123031e-02
## 119 Gammaproteobacteria South  214 96   1.118578e-01 6.899919e-04
## 120 Gammaproteobacteria Point  185 96   2.358815e-02 6.854785e-03
## 121 Gammaproteobacteria North  206 96   1.401483e-01 8.796677e-02
## 122 Gammaproteobacteria South  185 96   1.473569e-01 1.267272e-02
## 123 Gammaproteobacteria North  214 96   5.338684e-02 2.197811e-02
## 124 Gammaproteobacteria South  199 96   1.415582e-01 1.291130e-02
## 125 Gammaproteobacteria Point  172 96   1.868465e-02 5.322918e-03
## 126 Gammaproteobacteria South  206 64   8.356680e-02 1.164243e-03
## 127    Gemmatimonadetes North  185  9   2.371800e-03 1.259232e-04
## 128    Gemmatimonadetes North  178  9   1.149633e-03 6.022823e-04
## 129    Gemmatimonadetes North  199  9   3.603814e-03 1.021487e-03
## 130    Gemmatimonadetes South  178  9   1.604865e-03 5.802559e-04
## 131    Gemmatimonadetes North  172  9   1.654567e-03 3.606342e-04
## 132    Gemmatimonadetes Point  206  6   1.161221e-03 2.510368e-04
## 133    Gemmatimonadetes South  172  9   3.806662e-04 4.711227e-05
## 134    Gemmatimonadetes Point  214  9   2.243058e-03 2.879720e-04
## 135    Gemmatimonadetes North  206  9   2.625621e-03 1.223114e-03
## 136    Gemmatimonadetes South  185  9   1.962271e-03 6.206918e-04
## 137    Gemmatimonadetes Point  199  9   3.175512e-03 3.025115e-04
## 138    Gemmatimonadetes Point  185  9   1.573836e-03 7.500936e-04
## 139    Gemmatimonadetes North  214  9   2.273586e-03 6.391758e-04
## 140    Gemmatimonadetes Point  178  9   3.384843e-04 1.196620e-04
## 141    Gemmatimonadetes South  214  9   4.509576e-03 5.411028e-04
## 142    Gemmatimonadetes Point  172  9   5.768216e-04 3.820511e-04
## 143    Gemmatimonadetes South  199  9   1.630632e-03 3.054107e-04
## 144    Gemmatimonadetes South  206  6   2.047974e-03 1.041129e-04
## 145          Nitrospira North  206  9   5.714958e-04 5.987598e-04
## 146          Nitrospira North  178  9   3.878024e-04 4.690600e-04
## 147          Nitrospira South  185  9   6.402792e-04 2.130564e-04
## 148          Nitrospira North  199  9   2.357900e-04 1.075190e-04
## 149          Nitrospira North  214  9   1.677540e-04 4.710066e-05
## 150          Nitrospira Point  214  9   1.616646e-04 2.537888e-05
## 151          Nitrospira North  172  9   3.169511e-04 2.814620e-04
## 152          Nitrospira North  185  9   5.609904e-04 2.084704e-04
## 153          Nitrospira South  178  9   2.540769e-04 1.072225e-04
## 154          Nitrospira South  199  9   4.051250e-04 1.059307e-04
## 155          Nitrospira Point  199  9   1.008245e-04 8.197479e-06
## 156          Nitrospira Point  185  9   6.389268e-05 1.088997e-05
## 157          Nitrospira Point  206  6   1.346975e-04 1.000739e-05
## 158          Nitrospira South  172  9   9.705593e-05 2.440956e-06
## 159          Nitrospira Point  172  9   3.432994e-05 2.594432e-05
## 160          Nitrospira South  206  6   7.387673e-04 2.688982e-04
## 161          Nitrospira Point  178  9   2.270542e-05 2.247835e-05
## 162          Nitrospira South  214  9   1.609872e-03 7.345991e-04
## 163      Planctomycetia North  206 15   6.284917e-03 3.417952e-03
## 164      Planctomycetia South  185 15   6.014139e-03 4.269422e-04
## 165      Planctomycetia North  178 15   1.197726e-02 2.345113e-03
## 166      Planctomycetia Point  214 15   4.633208e-03 2.997879e-04
## 167      Planctomycetia North  214 15   4.458024e-03 3.959592e-04
## 168      Planctomycetia South  199 15   2.440070e-03 2.442694e-04
## 169      Planctomycetia Point  199 15   7.136162e-03 1.254046e-03
## 170      Planctomycetia North  199 15   3.438498e-03 8.727654e-04
## 171      Planctomycetia North  185 15   1.366957e-02 3.649779e-04
## 172      Planctomycetia South  172 15   2.979045e-03 3.465324e-04
## 173      Planctomycetia North  172 15   1.273890e-02 5.494993e-03
## 174      Planctomycetia Point  172 15   5.133297e-03 2.856916e-03
## 175      Planctomycetia South  206 10   8.466059e-03 1.415956e-03
## 176      Planctomycetia Point  185 15   1.111016e-02 3.365407e-03
## 177      Planctomycetia Point  178 15   5.062253e-03 1.037653e-03
## 178      Planctomycetia South  178 15   8.926841e-03 3.361142e-03
## 179      Planctomycetia Point  206 10   4.529090e-03 5.896583e-04
## 180      Planctomycetia South  214 15   4.943619e-03 2.196765e-04
## 181    Sphingobacteriia North  206  9   1.362852e-02 1.022626e-02
## 182    Sphingobacteriia North  178  9   4.609560e-03 1.324738e-03
## 183    Sphingobacteriia North  214  9   6.308727e-03 3.756637e-04
## 184    Sphingobacteriia South  199  9   1.483583e-02 7.841836e-04
## 185    Sphingobacteriia South  185  9   1.011084e-02 9.706777e-04
## 186    Sphingobacteriia North  185  9   4.101797e-03 6.293195e-04
## 187    Sphingobacteriia South  172  9   3.727195e-03 5.446491e-04
## 188    Sphingobacteriia Point  214  9   5.374118e-03 3.115481e-04
## 189    Sphingobacteriia Point  172  9   4.794169e-03 1.751406e-03
## 190    Sphingobacteriia South  206  6   7.213302e-03 1.240032e-03
## 191    Sphingobacteriia Point  199  9   5.938244e-03 1.142577e-03
## 192    Sphingobacteriia Point  185  9   4.027671e-03 9.521056e-04
## 193    Sphingobacteriia North  199  9   1.056244e-02 5.766108e-04
## 194    Sphingobacteriia South  178  9   4.200576e-03 1.372722e-03
## 195    Sphingobacteriia North  172  9   7.461864e-03 2.083922e-03
## 196    Sphingobacteriia Point  206  6   5.844065e-03 5.230740e-04
## 197    Sphingobacteriia Point  178  9   5.008239e-03 4.227298e-04
## 198    Sphingobacteriia South  214  9   1.409589e-02 2.554491e-03
## 199    Verrucomicrobiae North  214  3   1.466408e-03 1.119816e-04
## 200    Verrucomicrobiae North  206  3   2.642318e-03 9.759171e-04
## 201    Verrucomicrobiae North  185  3   1.033660e-03 2.257070e-04
## 202    Verrucomicrobiae North  178  3   1.689420e-03 3.296375e-04
## 203    Verrucomicrobiae Point  172  3   1.316726e-03 5.433870e-04
## 204    Verrucomicrobiae South  206  2   5.189818e-03 3.746604e-04
## 205    Verrucomicrobiae South  199  3   7.199899e-03 5.185931e-04
## 206    Verrucomicrobiae South  185  3   2.216577e-03 4.625670e-04
## 207    Verrucomicrobiae North  199  3   2.489547e-03 3.027166e-04
## 208    Verrucomicrobiae South  178  3   1.688051e-03 6.801064e-04
## 209    Verrucomicrobiae South  172  3   1.572256e-03 6.397107e-04
## 210    Verrucomicrobiae Point  214  3   3.397857e-03 3.632050e-04
## 211    Verrucomicrobiae North  172  3   2.973987e-03 5.876502e-04
## 212    Verrucomicrobiae Point  206  2   3.617375e-03 4.305849e-04
## 213    Verrucomicrobiae Point  199  3   2.771544e-03 3.021012e-04
## 214    Verrucomicrobiae Point  185  3   1.713061e-03 2.346569e-04
## 215    Verrucomicrobiae Point  178  3   8.660193e-04 2.424089e-04
## 216    Verrucomicrobiae South  214  3   5.493034e-03 2.367490e-04
##               se           ci ClassAbundance_zulkifly
## 1   2.953698e-03 5.894013e-03                    0.17
## 2   9.096768e-04 1.832183e-03                    0.17
## 3   1.307159e-03 2.608396e-03                    0.17
## 4   9.877412e-04 1.971007e-03                    0.17
## 5   6.461465e-04 1.289365e-03                    0.17
## 6   9.386353e-04 1.873018e-03                    0.17
## 7   9.452552e-04 1.886227e-03                    0.17
## 8   8.578237e-04 1.711760e-03                    0.17
## 9   2.939549e-03 5.865778e-03                    0.17
## 10  5.677182e-04 1.132864e-03                    0.17
## 11  1.451380e-03 2.896184e-03                    0.17
## 12  3.172443e-04 6.330512e-04                    0.17
## 13  1.565117e-03 3.123142e-03                    0.17
## 14  1.074124e-04 2.143381e-04                    0.17
## 15  9.204748e-04 1.853931e-03                    0.17
## 16  6.193033e-04 1.235801e-03                    0.17
## 17  1.501434e-03 2.996065e-03                    0.17
## 18  3.757978e-04 7.498929e-04                    0.17
## 19  7.060695e-04 1.451347e-03                    0.01
## 20  5.289170e-05 1.115917e-04                    0.01
## 21  3.741693e-05 7.691160e-05                    0.01
## 22  5.129623e-05 1.054409e-04                    0.01
## 23  5.487829e-05 1.128039e-04                    0.01
## 24  2.043322e-05 4.200109e-05                    0.01
## 25  2.482708e-05 5.103279e-05                    0.01
## 26  1.824720e-04 3.750766e-04                    0.01
## 27  3.016165e-05 6.199815e-05                    0.01
## 28  3.280636e-05 6.743443e-05                    0.01
## 29  1.613009e-05 3.315588e-05                    0.01
## 30  4.089616e-05 8.406327e-05                    0.01
## 31  3.144225e-05 6.463046e-05                    0.01
## 32  3.047418e-05 6.264058e-05                    0.01
## 33  3.070322e-05 6.311137e-05                    0.01
## 34  6.173464e-05 1.268974e-04                    0.01
## 35  3.469363e-05 7.319717e-05                    0.01
## 36  5.346739e-05 1.099038e-04                    0.01
## 37  1.886618e-03 3.810103e-03                    0.21
## 38  2.364816e-03 4.775843e-03                    0.21
## 39  3.247782e-03 6.559029e-03                    0.21
## 40  8.156239e-04 1.673522e-03                    0.21
## 41  1.434060e-03 2.896144e-03                    0.21
## 42  5.123066e-03 1.034624e-02                    0.21
## 43  5.563827e-04 1.123638e-03                    0.21
## 44  5.444398e-04 1.099518e-03                    0.21
## 45  7.259252e-04 1.466036e-03                    0.21
## 46  9.269902e-04 1.872095e-03                    0.21
## 47  7.507561e-04 1.516183e-03                    0.21
## 48  1.955722e-03 3.949661e-03                    0.21
## 49  1.847305e-03 3.730708e-03                    0.21
## 50  3.441187e-03 6.949619e-03                    0.21
## 51  2.643260e-03 5.338171e-03                    0.21
## 52  1.441657e-03 2.958036e-03                    0.21
## 53  1.247247e-03 2.518866e-03                    0.21
## 54  4.262855e-03 8.609011e-03                    0.21
## 55  4.130203e-03 9.524266e-03                    0.03
## 56  3.305931e-04 8.498167e-04                    0.03
## 57  4.474445e-03 1.031809e-02                    0.03
## 58  2.982661e-04 6.878027e-04                    0.03
## 59  5.161009e-04 1.190131e-03                    0.03
## 60  7.562419e-04 1.743897e-03                    0.03
## 61  1.494388e-03 3.446065e-03                    0.03
## 62  1.310787e-03 3.022680e-03                    0.03
## 63  4.026484e-03 9.285089e-03                    0.03
## 64  6.316763e-03 1.456648e-02                    0.03
## 65  1.838613e-03 4.239849e-03                    0.03
## 66  6.380940e-04 1.471447e-03                    0.03
## 67  5.897162e-03 1.359888e-02                    0.03
## 68  3.558353e-03 8.205577e-03                    0.03
## 69  2.975109e-04 6.860614e-04                    0.03
## 70  3.594133e-03 8.288085e-03                    0.03
## 71  4.930876e-05 1.267522e-04                    0.03
## 72  3.952121e-03 9.113608e-03                    0.03
## 73  1.329557e-04 2.647487e-04                    0.07
## 74  1.433305e-04 2.854075e-04                    0.07
## 75  1.126009e-04 2.242170e-04                    0.07
## 76  1.370534e-04 2.729082e-04                    0.07
## 77  2.538193e-04 5.095636e-04                    0.07
## 78  2.298777e-04 4.577451e-04                    0.07
## 79  1.614461e-04 3.214802e-04                    0.07
## 80  1.278769e-04 2.546355e-04                    0.07
## 81  5.044755e-04 1.004539e-03                    0.07
## 82  2.083215e-04 4.148212e-04                    0.07
## 83  3.216333e-04 6.404538e-04                    0.07
## 84  1.421383e-04 2.830335e-04                    0.07
## 85  1.886576e-04 3.756652e-04                    0.07
## 86  1.629886e-04 3.245517e-04                    0.07
## 87  1.867877e-03 3.719419e-03                    0.07
## 88  9.099977e-05 1.812037e-04                    0.07
## 89  1.094356e-04 2.197012e-04                    0.07
## 90  2.395453e-04 4.769957e-04                    0.07
## 91  8.875555e-04 1.953497e-03                    0.07
## 92  7.425314e-04 1.634301e-03                    0.07
## 93  2.042841e-03 4.496263e-03                    0.07
## 94  4.047494e-03 8.908475e-03                    0.07
## 95  1.267933e-03 2.790701e-03                    0.07
## 96  1.122077e-03 2.469675e-03                    0.07
## 97  1.047821e-03 2.306238e-03                    0.07
## 98  9.323323e-03 2.052050e-02                    0.07
## 99  7.790829e-04 1.714750e-03                    0.07
## 100 1.774728e-03 4.196565e-03                    0.07
## 101 1.219370e-03 2.683816e-03                    0.07
## 102 1.395915e-03 3.072388e-03                    0.07
## 103 7.702531e-04 1.695316e-03                    0.07
## 104 1.408848e-03 3.100853e-03                    0.07
## 105 8.884365e-04 1.955436e-03                    0.07
## 106 1.323509e-04 2.913024e-04                    0.07
## 107 5.469144e-04 1.203750e-03                    0.07
## 108 5.186076e-04 1.226312e-03                    0.07
## 109 3.279579e-04 6.510788e-04                    0.12
## 110 4.594023e-05 9.120289e-05                    0.12
## 111 3.520553e-05 6.989181e-05                    0.12
## 112 1.540345e-04 3.057971e-04                    0.12
## 113 2.275297e-04 4.517036e-04                    0.12
## 114 2.416966e-03 4.798284e-03                    0.12
## 115 2.180076e-03 4.327999e-03                    0.12
## 116 2.377835e-04 4.720599e-04                    0.12
## 117 6.695565e-04 1.329238e-03                    0.12
## 118 1.403789e-03 2.805248e-03                    0.12
## 119 7.042201e-05 1.398054e-04                    0.12
## 120 6.996135e-04 1.388908e-03                    0.12
## 121 8.978070e-03 1.782372e-02                    0.12
## 122 1.293404e-03 2.567731e-03                    0.12
## 123 2.243132e-03 4.453180e-03                    0.12
## 124 1.317754e-03 2.616072e-03                    0.12
## 125 5.432681e-04 1.078524e-03                    0.12
## 126 1.455304e-04 2.908192e-04                    0.12
## 127 4.197441e-05 9.679315e-05                    0.01
## 128 2.007608e-04 4.629551e-04                    0.01
## 129 3.404956e-04 7.851842e-04                    0.01
## 130 1.934186e-04 4.460242e-04                    0.01
## 131 1.202114e-04 2.772080e-04                    0.01
## 132 1.024854e-04 2.634470e-04                    0.01
## 133 1.570409e-05 3.621370e-05                    0.01
## 134 9.599067e-05 2.213549e-04                    0.01
## 135 4.077046e-04 9.401686e-04                    0.01
## 136 2.068973e-04 4.771060e-04                    0.01
## 137 1.008372e-04 2.325309e-04                    0.01
## 138 2.500312e-04 5.765729e-04                    0.01
## 139 2.130586e-04 4.913140e-04                    0.01
## 140 3.988734e-05 9.198036e-05                    0.01
## 141 1.803676e-04 4.159284e-04                    0.01
## 142 1.273504e-04 2.936705e-04                    0.01
## 143 1.018036e-04 2.347595e-04                    0.01
## 144 4.250390e-05 1.092598e-04                    0.01
## 145 1.995866e-04 4.602475e-04                    0.01
## 146 1.563533e-04 3.605514e-04                    0.01
## 147 7.101880e-05 1.637696e-04                    0.01
## 148 3.583965e-05 8.264639e-05                    0.01
## 149 1.570022e-05 3.620477e-05                    0.01
## 150 8.459628e-06 1.950794e-05                    0.01
## 151 9.382065e-05 2.163508e-04                    0.01
## 152 6.949014e-05 1.602445e-04                    0.01
## 153 3.574082e-05 8.241847e-05                    0.01
## 154 3.531025e-05 8.142558e-05                    0.01
## 155 2.732493e-06 6.301140e-06                    0.01
## 156 3.629989e-06 8.370770e-06                    0.01
## 157 4.085500e-06 1.050211e-05                    0.01
## 158 8.136520e-07 1.876285e-06                    0.01
## 159 8.648108e-06 1.994257e-05                    0.01
## 160 1.097772e-04 2.821913e-04                    0.01
## 161 7.492785e-06 1.727839e-05                    0.01
## 162 2.448664e-04 5.646628e-04                    0.01
## 163 8.825115e-04 1.892799e-03                    0.05
## 164 1.102360e-04 2.364327e-04                    0.05
## 165 6.055055e-04 1.298680e-03                    0.05
## 166 7.740489e-05 1.660170e-04                    0.05
## 167 1.022362e-04 2.192749e-04                    0.05
## 168 6.307009e-05 1.352719e-04                    0.05
## 169 3.237934e-04 6.944677e-04                    0.05
## 170 2.253470e-04 4.833213e-04                    0.05
## 171 9.423689e-05 2.021180e-04                    0.05
## 172 8.947428e-05 1.919032e-04                    0.05
## 173 1.418801e-03 3.043026e-03                    0.05
## 174 7.376524e-04 1.582107e-03                    0.05
## 175 4.477645e-04 1.012914e-03                    0.05
## 176 8.689445e-04 1.863701e-03                    0.05
## 177 2.679210e-04 5.746333e-04                    0.05
## 178 8.678432e-04 1.861339e-03                    0.05
## 179 1.864663e-04 4.218161e-04                    0.05
## 180 5.672022e-05 1.216528e-04                    0.05
## 181 3.408754e-03 7.860601e-03                    0.16
## 182 4.415793e-04 1.018284e-03                    0.16
## 183 1.252212e-04 2.887607e-04                    0.16
## 184 2.613945e-04 6.027769e-04                    0.16
## 185 3.235592e-04 7.461289e-04                    0.16
## 186 2.097732e-04 4.837378e-04                    0.16
## 187 1.815497e-04 4.186543e-04                    0.16
## 188 1.038494e-04 2.394771e-04                    0.16
## 189 5.838021e-04 1.346250e-03                    0.16
## 190 5.062408e-04 1.301333e-03                    0.16
## 191 3.808590e-04 8.782624e-04                    0.16
## 192 3.173685e-04 7.318532e-04                    0.16
## 193 1.922036e-04 4.432223e-04                    0.16
## 194 4.575741e-04 1.055168e-03                    0.16
## 195 6.946408e-04 1.601845e-03                    0.16
## 196 2.135441e-04 5.489325e-04                    0.16
## 197 1.409099e-04 3.249389e-04                    0.16
## 198 8.514971e-04 1.963556e-03                    0.16
## 199 6.465261e-05 2.781777e-04                    0.04
## 200 5.634460e-04 2.424312e-03                    0.04
## 201 1.303120e-04 5.606872e-04                    0.04
## 202 1.903163e-04 8.188650e-04                    0.04
## 203 3.137247e-04 1.349848e-03                    0.04
## 204 2.649249e-04 3.366191e-03                    0.04
## 205 2.994099e-04 1.288257e-03                    0.04
## 206 2.670632e-04 1.149080e-03                    0.04
## 207 1.747735e-04 7.519898e-04                    0.04
## 208 3.926596e-04 1.689478e-03                    0.04
## 209 3.693371e-04 1.589129e-03                    0.04
## 210 2.096965e-04 9.022513e-04                    0.04
## 211 3.392800e-04 1.459804e-03                    0.04
## 212 3.044695e-04 3.868652e-03                    0.04
## 213 1.744182e-04 7.504609e-04                    0.04
## 214 1.354792e-04 5.829200e-04                    0.04
## 215 1.399549e-04 6.021772e-04                    0.04
## 216 1.366871e-04 5.881170e-04                    0.04

p <- ggplot(relabund.classzulk.class.est.merge, aes(x=Date, y=ClassAbundance, color = Site)) + geom_point() +  geom_errorbar(aes(ymin=ClassAbundance-se, ymax=ClassAbundance+se)) +  geom_line() + facet_wrap(~Class, ncol = 3, scales="free_y") + scale_colour_hue(h=c(400, 120)) + geom_hline(aes(yintercept = ClassAbundance_zulkifly), linetype="dashed")+ scale_colour_hue(h=c(400, 120))
p + theme_bw() + theme(axis.text.x = element_text(size = 10, angle = 45, hjust=1),axis.text.y = element_text(size = 10))

# Summary of class abundance of classes_zulk of interest. 
relabund.classzulk.class.est.nosplit <- summarySE(relabund.classzulk.class, measurevar= "ClassAbundance", groupvars=c("Class"))
head(relabund.classzulk.class.est.nosplit)

##                 Class    N ClassAbundance          sd           se
## 1 Alphaproteobacteria 1196    0.118984942 0.026968565 7.798162e-04
## 2         Bacteroidia  468    0.001604779 0.001619960 7.488269e-05
## 3  Betaproteobacteria  728    0.139416618 0.032901155 1.219398e-03
## 4          Deinococci  156    0.059034453 0.023901398 1.913643e-03
## 5 Deltaproteobacteria 1352    0.019910446 0.009912747 2.695912e-04
## 6      Flavobacteriia  208    0.049475842 0.021336332 1.479408e-03
##             ci
## 1 0.0015299613
## 2 0.0001471487
## 3 0.0023939616
## 4 0.0037801865
## 5 0.0005288628
## 6 0.0029166396

relabund.classzulk.class.est.merge.nosplit <- merge(relabund.classzulk.class.est.nosplit, classzulk, by ="Class")
relabund.classzulk.class.est.merge.nosplit

##                  Class    N ClassAbundance           sd           se
## 1  Alphaproteobacteria 1196   0.1189849423 0.0269685654 7.798162e-04
## 2          Bacteroidia  468   0.0016047794 0.0016199603 7.488269e-05
## 3   Betaproteobacteria  728   0.1394166181 0.0329011553 1.219398e-03
## 4           Deinococci  156   0.0590344529 0.0239013981 1.913643e-03
## 5  Deltaproteobacteria 1352   0.0199104462 0.0099127474 2.695912e-04
## 6       Flavobacteriia  208   0.0494758417 0.0213363322 1.479408e-03
## 7  Gammaproteobacteria 1664   0.0668946728 0.0535523210 1.312809e-03
## 8     Gemmatimonadetes  156   0.0019508196 0.0012040756 9.640320e-05
## 9           Nitrospira  156   0.0003584377 0.0004587545 3.672975e-05
## 10      Planctomycetia  260   0.0069005428 0.0040347804 2.502265e-04
## 11    Sphingobacteriia  156   0.0073552263 0.0044809278 3.587613e-04
## 12    Verrucomicrobiae   52   0.0026770285 0.0017143254 2.377342e-04
##              ci ClassAbundance_zulkifly
## 1  1.529961e-03                    0.17
## 2  1.471487e-04                    0.01
## 3  2.393962e-03                    0.21
## 4  3.780187e-03                    0.03
## 5  5.288628e-04                    0.07
## 6  2.916640e-03                    0.07
## 7  2.574933e-03                    0.12
## 8  1.904336e-04                    0.01
## 9  7.255548e-05                    0.01
## 10 4.927373e-04                    0.05
## 11 7.086925e-04                    0.16
## 12 4.772712e-04                    0.04

df <- relabund.classzulk.class.est.merge.nosplit
p <- ggplot(df, aes(Class)) 
p + geom_point(aes(y = df$ClassAbundance, color = "Braus 2016 (2014)")) + 
  geom_point(aes(y = df$ClassAbundance_zulkifly, color = "Zulkifly 2012 (2011)")) + 
  labs(title ="Cladophora, Lake Mendota") +
  ylab("Relative Abundance")+ 
  guides(color=guide_legend(title="Growth Season")) + 
  theme_bw() + 
  theme(axis.text.x = element_text(size = 10, angle = 55, hjust=1),axis.text.y = element_text(size = 10))